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Abstract 

Gradient-based approaches to direct policy search in reinforcement learning have received 
much recent attention as a means to solve problems of partial observability and to avoid some of 
the problems associated with policy degradation in value-function methods. In this paper we intro- 
duce GPOMDP, a simulation-based algorithm for generating a biased estimate of the gradient of 
the average reward in Partially Observable Markov Decision Processes (POMDPs) controlled by 
parameterized stochastic policies. A similar algorithm was proposed by Kimura, Yamamura, and 
Kobayashi (1995). The algorithm's chief advantages are that it requires storage of only twice the 
number of policy parameters, uses one free parameter /3 G [0, 1) (which has a natural interpretation 
in terms of bias-variance trade-off), and requires no knowledge of the underlying state. We prove 
convergence of GPOMDP, and show how the correct choice of the parameter /? is related to the 
mixing time of the controlled POMDP. We briefly describe extensions of GPOMDP to controlled 
Markov chains, continuous state, observation and control spaces, multiple-agents, higher-order 
derivatives, and a version for training stochastic policies with internal states. In a companion paper 
(Baxter, Bartlett, & Weaver, 2001) we show how the gradient estimates generated by GPOMDP 
can be used in both a traditional stochastic gradient algorithm and a conjugate-gradient procedure 
to find local optima of the average reward. 

1. Introduction 

Dynamic Programming is the method of choice for solving problems of decision making under 
uncertainty (Bertsekas, 1995). However, the application of Dynamic Programming becomes prob- 
lematic in large or infinite state-spaces, in situations where the system dynamics are unknown, or 
when the state is only partially observed. In such cases one looks for approximate techniques that 
rely on simulation, rather than an explicit model, and parametric representations of either the value- 
function or the policy, rather than exact representations. 

Simulation-based methods that rely on a parametric form of the value function tend to go by 
the name "Reinforcement Learning," and have been extensively studied in the Machine Learning 
literature (Bertsekas & Tsitsiklis, 1996; Sutton & Barto, 1998). This approach has yielded some 
remarkable empirical successes in a number of different domains, including learning to play check- 
ers (Samuel, 1959), backgammon (Tesauro, 1992, 1994), and chess (Baxter, Tridgell, & Weaver, 
2000), job-shop scheduling (Zhang & Dietterich, 1995) and dynamic channel allocation (Singh & 
Bertsekas, 1997). 

Despite this success, most algorithms for training approximate value functions suffer from the 
same theoretical flaw: the performance of the greedy policy derived from the approximate value- 
function is not guaranteed to improve on each iteration, and in fact can be worse than the old policy 
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by an amount equal to the maximum approximation error over all states. This can happen even when 
the parametric class contains a value function whose corresponding greedy policy is optimal. We 
illustrate this with a concrete and very simple example in Appendix A. 

An alternative approach that circumvents this problem — the approach we pursue here — is to 
consider a class of stochastic policies parameterized by € M. K , compute the gradient with respect 
to of the average reward, and then improve the policy by adjusting the parameters in the gradient 
direction. Note that the policy could be directly parameterized, or it could be generated indirectly 
from a value function. In the latter case the value-function parameters are the parameters of the 
policy, but instead of being adjusted to minimize error between the approximate and true value 
function, the parameters are adjusted to directly improve the performance of the policy generated 
by the value function. 

These "policy-gradient" algorithms have a long history in Operations Research, Statistics, Con- 
trol Theory, Discrete Event Systems and Machine Learning. Before describing the contribution of 
the present paper, it seems appropriate to introduce some background material explaining this ap- 
proach. Readers already familial - with this material may want to skip directly to section 1.2, where 
the contributions of the present paper are described. 

1.1 A Brief History of Policy-Gradient Algorithms 

For large-scale problems or problems where the system dynamics are unknown, the performance 
gradient will not be computable in closed form 1 . Thus the challenging aspect of the policy-gradient 
approach is to find an algorithm for estimating the gradient via simulation. Naively, the gradient 
can be calculated numerically by adjusting each parameter in turn and estimating the effect on per- 
formance via simulation (the so-called crude Monte-Carlo technique), but that will be prohibitively 
inefficient for most problems. Somewhat surprisingly, under mild regularity conditions, it turns out 
that the full gradient can be estimated from a single simulation of the system. The technique is 
called the score function or likelihood ratio method and appears to have been first proposed in the 
sixties (Aleksandrov, Sysoyev, & Shemeneva, 1968; Rubinstein, 1969) for computing performance 
gradients in i.i.d. (independently and identically distributed) processes. 

Specifically, suppose r(X) is a performance function that depends on some random variable 
X, and q(0,x) is the probability that X = x, parameterized by G K . Under mild regularity 
conditions, the gradient with respect to of the expected performance, 



T]{6) = Er(X) 



(1) 



may be written 



Er(X) 



Vg(MQ 
q(B,X) ' 



(2) 



To see this, rewrite (1) as a sum 



77(0) 



X 



differentiate (one source of the requirement of "mild regularity conditions") to obtain 



Vrj(0) = J2r(x)Vq(0,x) 



X 



1. See equation (17) for a closed-form expression for the performance gradient. 
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rewrite as 

and observe that this formula is equivalent to (2). 

If a simulator is available to generate samples X distributed according to q(8,x), then any 
sequence X±, X2, ■ ■ ■ , X^ generated i.i.d. according to q{6 : x) gives an unbiased estimate, 

of Vr}{6). By the law of large numbers, Vr/(#) — > Vr/(#) with probability one. The quantity 
Vq(6,X)/q(6,X) is known as the likelihood ratio or score function in classical statistics. If 
the performance function r(X) also depends on 6, then r(X)Vq(6,X)/q(6,X) is replaced by 
Vr(0, X) + r(0, X)Vq(0, X)/q{0, X) in (2). 

1.1.1 Unbiased Estimates of the Performance Gradient for Regenerative 
Processes 

Extensions of the likelihood-ratio method to regenerative processes (including Markov Decision 
Processes or MDPs) were given by Glynn (1986, 1990), Glynn and L'Ecuyer (1995) and Reiman 
and Weiss (1986, 1989), and independently for episodic Partially Observable Markov Decision 
Processes (POMDPs) by Williams (1992), who introduced the REINFORCE algorithm 2 . Here the 
i.i.d. samples X of the previous section are sequences of states Xq,... ,Xt (of random length) 
encountered between visits to some designated recurrent state i*, or sequences of states from some 
start state to a goal state. In this case Vq(<9, X)/q(6, X) can be written as a sum 

q{0,X) ^ Px t x t+ M ' 

where px t x t+ i (#) lS the transition probability from X t to X t+ \ given parameters 6. Equation (4) 
admits a recursive computation over the course of a regenerative cycle of the form z = £ M. K , 
and after each state transition X t — > X t+ \, 

, Vpx t x t+1 (0) _ 

z t+ i = z t H 7^—, (5) 

Px t x t+1 [0) 

so that each term r(X)Vq(0, X)/q(6, X) in the estimate (3) is of the form 3 t{Xq, . . . , Xt)zt- If, 
in addition, t(Xq, . . . , Xt) can be recursively computed by 

r{X , . . . , X t+1 ) = 4>{r{X , X t ),X t+1 ) 

for some function <fi, then the estimate r(Xo, ■ ■ ■ ,Xt)zt for each cycle can be computed using 
storage of only K + 1 parameters (K for z% and 1 parameter to update the performance function 
r). Hence, the entire estimate (3) can be computed with storage of only 2K + 1 real parameters, as 
follows. 



2. A thresholded version of these algorithms for neuron-like elements was described earlier in Barto, Sutton, and An- 
derson (1983). 

3. The vector zt is known in reinforcement learning as an eligibility trace. This terminology is used in Barto et al. 
(1983). 
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Algorithm 1.1: Policy-Gradient Algorithm for Regenerative Processes. 

1. Set j = 0, r = 0, z = 0, and A = (z , A G M K ). 

2. For each state transition X t — >■ Xt+i: 

• If the episode is finished (that is, X t+ \ = i*), set 
A j+1 = Aj + r t z t , 

3=3 + 1> 
= 0, 

n+i = 0. 

• Otherwise, set 

+ PXtXt+i(fl) , 
n+i = (f>(r t ,X t +i). 

3. If j = N return A^/N, otherwise goto 2. 

Examples of recursive performance functions include the sum of a scalar reward over a cycle, 
r(Xo, . . . ,Xt) = Y^d=o r {Xt) where r(i) is a scalar reward associated with state i (this corre- 
sponds to rj(8) being the average reward multiplied by the expected recurrence time E# [T]); the 
negative length of the cycle (which can be implemented by assigning a reward of — 1 to each state, 
and is used when the task is to mimimize time taken to get to a goal state, since r](6) in this case is 
just —Eg [T]); the discounted reward from the start state, r(Xo, . . . , Xt) = YlJ=o «M^*)> where 
a 6 [0, 1) is the discount factor, and so on. 

As Williams (1992) pointed out, a further simplification is possible in the case that tt = 
r(Xo, . . . , Xt) is a sum of scalar rewards r(X t , t) depending on the state and possibly the time 
t since the starting state (such as r(Xt , t) = r(X t ), or r(Xt, t) = a t r{X t ) as above). In that case, 
the update A from a single regenerative cycle may be written as 



A = y^ 1 Vpx t x t+1 (0) 



t=0 



,s=0 s=t+l 



Because changes in px t Xt+i(@) have no influence on the rewards r(X s , s) associated with earlier 
states (s < t), we should be able to drop the first term in the parentheses on the right-hand-side and 
write 



A= L pxx W) L r ( X »*)- < 6 > 
t=0 PXtXt+iK ) s=t+ i 



Although the proof is not entirely trivial, this intuition can indeed be shown to be correct. 

Equation (6) allows an even simpler recursive formula for estimating the performance gradi- 
ent. Set zo = Aq = 0, and introduce a new variable s = 0. As before, set zt+i = Zt + 
'Vpx t Xt + i(8)/Px t Xt +1 (8) an d s = s + 1 if X t +i / %*, or s = and z t +i = otherwise. But 
now, on each iteration, set A i+ i = r(X t , s)z t + A*. Then At/t is our estimate of Vr)(6). Since At 
is updated on every iteration, this suggests that we can do away with At altogether and simply up- 
date 6 directly: = #t+7t?"(^t, s)z t , where the j t are suitable step-sizes 4 . Proving convergence 

4. The usual requirements on 74 for convergence of a stochastic gradient algorithm are "ft > 0, X^ilo 7* = °°> an ^ 

Eoo 2 ^ 
4=0 It < °°- 
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of such an algorithm is not as straightforward as normal stochastic gradient algorithms because the 
updates r(X t )z t are not in the gradient direction (in expectation), although the sum of these updates 
over a regenerative cycle are. Marbach and Tsitsiklis (1998) provide the only convergence proof that 
we know of, albeit for a slightly different update of the form 6 t +\ = t + jt [ r {Xt, s) — f){Qt)) z u 
where fj(9t) is a moving estimate of the expected performance, and is also updated on-line (this 
update was first suggested in the context of POMDPs by Jaakkola et al. (1995)). 

Marbach and Tsitsiklis (1998) also considered the case of ^-dependent rewards (recall the dis- 
cussion after (3)), as did Baird and Moore (1999) with their "VAPS" algorithm (Value And Policy 
Search). This last paper contains an interesting insight: through suitable choices of the performance 
function r(Xo, ■ ■ ■ , Xt, 9), one can combine policy-gradient search with approximate value func- 
tion methods. The resulting algorithms can be viewed as actor-critic techniques in the spirit of Barto 
et al. (1983); the policy is the actor and the value function is the critic. The primary motivation is 
to reduce variance in the policy-gradient estimates. Experimental evidence for this phenomenon 
has been presented by a number of authors, including Barto et al. (1983), Kimura and Kobayashi 
(1998a), and Baird and Moore (1999). More recent work on this subject includes that of Sutton 
et al. (2000) and Konda and Tsitsiklis (2000). We discuss the use of VAPS-style updates further in 
Section 6.2. 

So far we have not addressed the question of how the parameterized state-transition probabili- 
ties pxtXt+i ($) arise. Of course, they could simply be generated by parameterizing the matrix of 
transition probabilities directly. Alternatively, in the case of MDPs or POMDPs, state transitions 
are typically generated by feeding an observation Y t that depends stochastically on the state X t 
into a parameterized stochastic policy, which selects a control U t at random from a set of avail- 
able controls (approximate value-function based approaches that generate controls stochastically 
via some form of lookahead also fall into this category). The distribution over successor states 
PXtXt+i {Ut) is then a fixed function of the control. If we denote the probability of control u t given 
parameters 9 and observation y t by fi Vt (9, yt), then all of the above discussion carries through with 
Vpx t Xt +1 (0)/pxtXt+i(d) replaced by V/j, Ut (9,Y t )/ij Ut (9,Yt). In that case, Algorithm 1.1 is pre- 
cisely Williams' REINFORCE algorithm. 

Algorithm 1.1 and the variants above have been extended to cover multiple agents (Peshkin 
et al., 2000), policies with internal state (Meuleau et al., 1999), and importance sampling methods 
(Meuleau et al., 2000). We also refer the reader to the work of Rubinstein and Shapiro (1993) 
and Rubinstein and Melamed (1998) for in-depth analysis of the application of the likelihood-ratio 
method to Discrete-Event Systems (DES), in particular networks of queues. Also worth mentioning 
is the large literature on Infinitesimal Perturbation Analysis (IPA), which seeks a similar goal of esti- 
mating performance gradients, but operates under more restrictive assumptions than the likelihood- 
ratio approach; see, for example, Ho and Cao (1991). 

1.1.2 Biased Estimates of the Performance Gradient 

All the algorithms described in the previous section rely on an identifiable recurrent state i*, either 
to update the gradient estimate, or in the case of the on-line algorithm, to zero the eligibility trace 
z. This reliance on a recurrent state can be problematic for two main reasons: 

1. The variance of the algorithms is related to the recurrence time between visits to i* , which 
will typically grow as the state space grows. Furthermore, the time between visits depends on 
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the parameters of the policy, and states that are frequently visited for the initial value of the 
parameters may become very rare as performance improves. 

2. In situations of partial observability it may be difficult to estimate the underlying states, and 
therefore to determine when the gradient estimate should be updated, or the eligibility trace 
zeroed. 

If the system is available only through simulation, it seems difficult (if not impossible) to obtain 
unbiased estimates of the gradient direction without access to a recurrent state. Thus, to solve 1 
and 2, we must look to biased estimates. Two principle techniques for introducing bias have been 
proposed, both of which may be viewed as artificial truncations of the eligibility trace z. The first 
method takes as a starting point the formula 5 for the eligibility trace at time t: 

= Vpx s x B+1 (0) 

and simply truncates it at some (fixed, not random) number of terms n looking backwards (Glynn, 
1990; Rubinstein, 1991, 1992; Cao & Wan, 1998): 

The eligibility trace Zt(n) is then updated after each transition X t — > X t+ i by 

z t+1 (n )=z t (n + — — , (8) 

PX t X t+1 (0) PXt-nXt-n+i (°) 

and in the case of state-based rewards r(X t ), the estimated gradient direction after T steps is 

t—n 

Unless n exceeds the maximum recurrence time (which is infinite in an ergodic Markov chain), 
V n r](6) is a biased estimate of the gradient direction, although as n — > oo, the bias approaches zero. 
However the variance of V n r]{0) diverges in the limit of large n. This illustrates a natural trade-off 
in the selection of the parameter n: it should be large enough to ensure the bias is acceptable (the 
expectation of V n 'q{9) should at least be within 90° of the true gradient direction), but not so large 
that the variance is prohibitive. Experimental results by Cao and Wan (1998) illustrate nicely this 
bias/variance trade-off. 

One potential difficulty with this method is that the likelihood ratios Vpx s x, +1 {0) /px s x, +1 {0) 
must be remembered for the previous n time steps, requiring storage of Kn parameters. Thus, 
to obtain small bias, the memory may have to grow without bound. An alternative approach that 
requires a fixed amount of memory is to discount the eligibility trace, rather than truncating it: 

z t+l {f3) :=pz t (P) + VPXtXt+l< f \ (10) 
Px t x t+1 (9) 



5. For ease of exposition, we have kept the expression for z in terms of the likelihood ratios Vpx s x s+1 {6) /px B x e+1 {8) 
which rely on the availability of the underlying state X s . If X s is not available, Vpx s x s+1 (0) /px,x s+1 (0) should 
be replaced with V ' ji Us (0, Y B )/fJ,u, (0, Y B ). 
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where zq(/3) = and /3 € [0, 1) is a discount factor. In this case the estimated gradient direction 
after T steps is simply 



This is precisely the estimate we analyze in the present paper. A similar estimate with r(X t )z t (/3) 
replaced by (r(X t ) — b)z t (P) where b is a reward baseline was proposed by Kimura et al. (1995, 
1997) and for continuous control by Kimura and Kobayashi (1998b). In fact the use of (r(X t ) — b) 
in place of r(X t ) does not affect the expectation of the estimates of the algorithm (although judi- 
cious choice of the reward baseline b can reduce the variance of the estimates). While the algorithm 
presented by Kimura et al. (1995) provides estimates of the expectation under the stationary distri- 
bution of the gradient of the discounted reward, we will show that these are in fact biased estimates 
of the gradient of the expected discounted reward. This arises because the stationary distribution 
itself depends on the parameters. A similar estimate to (11) was also proposed by Marbach and 
Tsitsiklis (1998), but this time with r(X t )z t (/3) replaced by (r(X t ) — fi(6))z t (/3), where fj(6) is an 
estimate of the average reward, and with zt zeroed on visits to an identifiable recurrent state. 

As a final note, observe that the eligibility traces z t {0) and z t (n) defined by (10) and (8) are 
simply filtered versions of the sequence Vpx t x t+1 (6)/px t x t +i ($)> a first-order, infinite impulse 
response filter in the case of z t (/3) and an ra-th order, finite impulse response filter in the case of 
zt(n). This raises the question, not addressed in this paper, of whether there is an interesting theory 
of optimal filtering for policy-gradient estimators. 

1.2 Our Contribution 

We describe GPOMDP, a general algorithm based upon (11) for generating a biased estimate of the 
performance gradient Vrj(6) in general POMDPs controlled by parameterized stochastic policies. 
Here r)(0) denotes the average reward of the policy with parameters 6 £ M. K . GPOMDP does 
not rely on access to an underlying recurrent state. Writing Vprj{6) for the expectation of the esti- 
mate produced by GPOMDP, we show that lim^! VpT]{6) = Vrj(6), and more quantitatively that 
VbT)(6) is close to the true gradient provided 1/(1-/3) exceeds the mixing time of the Markov chain 
induced by the POMDP 6 . As with the truncated estimate above, the trade-off preventing the setting 
of j3 arbitrarily close to 1 is that the variance of the algorithm's estimates increase as f} approaches 
1. We prove convergence with probability 1 of GPOMDP for both discrete and continuous observa- 
tion and control spaces. We present algorithms for both general parameterized Markov chains and 
POMDPs controlled by parameterized stochastic policies. 

There are several extensions to GPOMDP that we have investigated since the first version of 
this paper was written. We outline these developments briefly in Section 7. 

In a companion paper we show how the gradient estimates produced by GPOMDP can be used 
to perform gradient ascent on the average reward r](8) (Baxter et al., 2001). We describe both 
traditional stochastic gradient algorithms, and a conjugate-gradient algorithm that utilizes gradient 
estimates in a novel way to perform line searches. Experimental results are presented illustrat- 

6. The mixing-time result in this paper applies only to Markov chains with distinct eigenvalues. Better estimates of the 
bias and variance of GPOMDP may be found in Bartlett and Baxter (2001), for more general Markov chains than 
those treated here, and for more refined notions of the mixing time. Roughly speaking, the variance of GPOMDP 
grows with 1/(1 — /?), while the bias decreases as a function of 1/(1 — j3). 




T-1 



(ID 



t=0 



325 



Baxter & Bartlett 



ing both the theoretical results of the present paper on a toy problem, and practical aspects of the 
algorithms on a number of more realistic problems. 

2. The Reinforcement Learning Problem 

We model reinforcement learning as a Markov decision process (MDP) with a finite state space 
S = {l,...,n}, and a stochastic matrix 7 P = [pjj] giving the probability of transition from state 
i to state j. Each state i has an associated reward 8 r(i). The matrix P belongs to a parameterized 
class of stochastic matrices, V := {P(8): 8 £ M. K }. Denote the Markov chain corresponding to 
P(9) by M{6). We assume that these Markov chains and rewards satisfy the following assumptions: 

Assumption 1. Each P{9) £ V has a unique stationary distribution tt(8) := [7r(#, 1), . . . , 7r(#, n)]' 
satisfying the balance equations 

n'(8)P(8) = ir'{0) (12) 
(throughout ir' denotes the transpose ofir). 

Assumption 2. The magnitudes of the rewards, \r(i)\, are uniformly bounded by R < oo for all 
states i. 

Assumption 1 ensures that the Markov chain forms a single recurrent class for all parameters 8. 
Since any finite-state Markov chain always ends up in a recurrent class, and it is the properties of 
this class that determine the long-term average reward, this assumption is mainly for convenience 
so that we do not have to include the recurrence class as a quantifier in our theorems. However, 
when we consider gradient-ascent algorithms Baxter et al. (2001), this assumption becomes more 
restrictive since it guarantees that the recurrence class cannot change as the parameters are adjusted. 

Ordinarily, a discussion of MDPs would not be complete without some mention of the actions 
available in each state and the space of policies available to the learner. In particular, the parameters 
8 would usually determine a policy (either directly or indirectly via a value function), which would 
then determine the transition probabilities P(8). However, for our puiposes we do not care how 
the dependence of P on 8 arises, just that it satisfies Assumption 1 (and some differentiability 
assumptions that we shall meet in the next section). Note also that it is easy to extend this setup 
to the case where the rewards also depend on the parameters 8 or on the transitions i — > j. It is 
equally straightforward to extend our algorithms and results to these cases. See Section 6.1 for an 
illustration. 

The goal is to find a 8 G M. K maximizing the average reward: 

X = i , 

where E# denotes the expectation over all sequences Xq , X\ , . . . , with transitions generated ac- 
cording to P(8). Under Assumption 1, rj(8) is independent of the starting state % and is equal to 

n 

v ( 6 ) = Y / 7r(8,i)r(i)=TT'(8)r, (13) 

2=1 

where r = [r(l), . . . , r(n)]' (Bertsekas, 1995). 

7. A stochastic matrix P = \pij] has pij > for all i, j and X^J=i Pij = 1 f° r a ^ *■ 

8. All the results in the present paper apply to bounded stochastic rewards, in which case r(i) is the expectation of the 
reward in state i. 



rj(8) := lim E 6 

T->oo 



T-l 
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3. Computing the Gradient of the Average Reward 



For general MDPs little will be known about the average reward rf(6), hence finding its optimum 
will be problematic. However, in this section we will see that under general assumptions the gradient 
Vr/(#) exists, and so local optimization of rj(8) is possible. 

To ensure the existence of suitable gradients (and the boundedness of certain random variables), 
we require that the parameterized class of stochastic matrices satisfies the following additional as- 
sumption. 



Assumption 3. The derivatives, 

VP(8) := 

exist for all 6 £ K . The ratios 



i,j=l...n;k=l...K 



d0 k 



are uniformly bounded by B < oo for all 9 6 



i,j—l...n;k—l...K 



The second part of this assumption allows zero-probability transitions pij(6) = only if 
Vpij{6) is also zero, in which case we set 0/0 = 0. One example is if i — > j is a forbidden 
transition, so that pjj (6) = for all 6 £ M. h . Another example satisfying the assumption is 



Pi 



e *■» 



E 



3=1 



where 0=[8 n ,.. .,0 ln , 



5 8nn] £ ar e the parameters of P{ff), for then 



P*j(0) 
d Pij (6)/d6 kl 

Pi 0) 



1 



-Pij(6), and 
PkiiO). 



Assuming for the moment that Vn(6) exists (this will be justified shortly), then, suppressing 9 
dependencies, 

Vr/ = Vvr'r, (14) 

since the reward r does not depend on 6. Note that our convention for V in this paper is that it takes 
precedence over all other operations, so Vgi0)fi9) = [Vg(#)j f(9). Equations like (14) should be 
regarded as shorthand notation for K equations of the form 



dr/(6>) 
00k 



0tt(0,1) d7r(0,n) 



[r(l),...,r(n)]' 



89 k ' ' 89 k 

where k = 1, . . . , K . To compute V7r, first differentiate the balance equations (12) to obtain 

Vtt'P + tt'VP = Vtt', 
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and hence 



V7r'(J - P) = tt'VP 



(15) 



The system of equations defined by (15) is under-constrained because / — P is not invertible (the 
balance equations show that I — P has a left eigenvector with zero eigenvalue). However, let e 
denote the n-dimensional column vector consisting of all Is, so that en' is the n x n matrix with the 
stationary distribution tt' in each row. Since W'e = V(7r'e) = V(l) = 0, we can rewrite (15) as 

Vtt' [I-(P- en')] = tt'VP 

To see that the inverse [I — (P — en')]^ 1 exists, let A be any matrix satisfying lim i _ ) . 00 A 1 = 0. 
Then we can write 



lim 

T-s-oo 



T 



t=o 



T 



T+l 



lim 

T-s-oo 

_t=0 t=l 

I- lim A T+1 



Thus, 



(i-a^^a*. 



t=0 



It is easy to prove by induction that [P — en'] 1 = P l - en' which converges to as t — > oo by 
Assumption 1. So [I — (P — en')] -1 exists and is equal to Yltlo ~ e7r '] • H ence > we can write 



Vtt' = tt'VP [i - P + en'] ' 



(16) 



and so 9 



Vr/ = tt'VP [i - P + en'] 1 r. 



(17) 



For MDPs with a sufficiently small number of states, (17) could be solved exactly to yield the precise 
gradient direction. However, in general, if the state space is small enough that an exact solution of 
(17) is possible, then it will be small enough to derive the optimal policy using policy iteration and 
table-lookup, and there would be no point in pursuing a gradient based approach in the first place 10 . 

Thus, for problems of practical interest, (17) will be intractable and we will need to find some 
other way of computing the gradient. One approximate technique for doing this is presented in the 
next section. 



9. The argument leading to (16) coupled with the fact that n(8) is the unique solution to (12) can be used to justify the 
existence of Vtt. Specifically, we can run through the same steps computing the value of ir(8 + S) for small d and 
show that the expression (16) for X7ir is the unique matrix satisfying tv(6 + 5) = n(8) + <5V-7r(#) + 0(||5|| 2 ). 
10. Equation (17) may still be useful for POMDPs, since in that case there is no tractable dynamic programming 
algorithm. 
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4. Approximating the Gradient in Parameterized Markov Chains 



In this section, we show that the gradient can be split into two components, one of which becomes 
negligible as a discount factor [3 approaches 1. 

For all /3 G [0, 1), let Jg(6) = [Jg(0, 1), . . . , Jg(0, n)] denote the vector of expected discounted 
rewards from each state i: 



Jg{0,i) :=E 



(18) 



Where the 6 dependence is obvious, we just write Jg. 
Proposition 1. For all G R K and j3 G [0, 1), 

Vr; = (1 - flVit'Jp + fin'VPJp. 
Proof. Observe that Jg satisfies the Bellman equations: 

Jg = r + pPJg. 

(Bertsekas, 1995). Hence, 

Vr/ = V7rV 

= Vtt' [J P - pPJ/s] 

= Vn'jg - fiVn'jg + fin'VPJp 

= (1 - P)Vn'jp + fin'VPJp. 



(19) 



(20) 



by (15) 



□ 



We shall see in the next section that the second term in (19) can be estimated from a single sam- 
ple path of the Markov chain. In fact, Theorem 1 in (Kimura et al., 1997) shows that the gradient 
estimates of the algorithm presented in that paper converge to (1 — /3)ir'V Jg. By the Bellman equa- 
tions (20), this is equal to (1 - /3)P(tt'VPJ b +tt'VJ b ), which implies (1 - P)tt'VJ 8 = Ptt'VPJb- 
Thus the algorithm of Kimura et al. (1997) also estimates the second term in the expression for 
Vr](6) given by (19). It is important to note that n'VJg / V [tt'Jb] — the two quantities disagree 
by the first term in (19). This arises because the the stationary distribution itself depends on the 
parameters. Hence, the algorithm of Kimura et al. (1997) does not estimate the gradient of the ex- 
pected discounted reward. In fact, the expected discounted reward is simply 1/(1 — 0) times the 
average reward rj(8) (Singh et al., 1994, Fact 7), so the gradient of the expected discounted reward 
is proportional to the gradient of the average reward. 

The following theorem shows that the first term in (19) becomes negligible as j3 approaches 1. 
Notice that this is not immediate from Proposition 1 , since Jg can become arbitrarily large in the 
limit 1. 



Theorem 2. For all 6 G 



where 



limV^, 



\7gr) := n'VPJg 



(21) 



(22) 
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Proof. Recalling equation (17) and the discussion preceeding it, we have 11 

oo 

V77 = tt'VP i pt - e7r '] r - 



(23) 



But VPe = V(Pe) = V(l) = since P is a stochastic matrix, so (23) can be rewritten as 



V?7 = 7r' 



]T VPP* 



t=o 



r. 



(24) 



Now let /3 E [0, 1] be a discount factor and consider the expression 



HP) ■■= 



t=0 



(25) 



Clearly Vrj = lim^i /(/3). To complete the proof we just need to show that = Vpr). 
Since (/3P)* = P t P t — > (3* en' — >■ 0, we can invoke the observation before (16) to write 



J2(pp) t = [i-ppy 



t=0 



In particular, 5^f. (/3P)* converges, so we can take VP back out of the sum in the right-hand-side 
of (25) and write 12 



HP) = ^'VP 



5>P* 



t=o 



But [JXo P'P 1 ] r = Jp. Thus f(P) = n'VPJ p = V^. 



(26) 
□ 



Theorem 2 shows that Vpr] is a good approximation to the gradient as f3 approaches 1, but it 
turns out that values of j3 very close to 1 lead to large variance in the estimates of \Z377 that we 
describe in the next section. However, the following theorem shows that 1 — (3 need not be too 
small, provided the transition probability matrix P(0) has distinct eigenvalues, and the Markov 
chain has a short mixing time. From any initial state, the distribution over states of a Markov chain 
converges to the stationary distribution, provided the assumption (Assumption 1) about the existence 
and uniqueness of the stationary distribution is satisfied (see, for example, Lancaster & Tismenetsky, 
1985, Theorem 15.8.1, p. 552). The spectral resolution theorem (Lancaster & Tismenetsky, 1985, 
Theorem 9.5.1, p. 314) implies that the distribution converges to stationarity at an exponential rate, 
and the time constant in this convergence rate (the mixing time) depends on the eigenvalues of 
the transition probability matrix. The existence of a unique stationary distribution implies that the 

11. Since e-n'r = erj, (23) motivates a different kind of algorithm for estimating V?? based on differential rewards 
(Marbach & Tsitsiklis, 1998). 

12. We cannot back VP out of the sum in the right-hand-side of (24) because diverges (P* — > e7r')- The reason 
^2^L VPP* converges is that P becomes orthogonal to VP in the limit of large t. Thus, we can view Y^tLo * 
as a sum of two orthogonal components: an infinite one in the direction e and a finite one in the direction e~ L . It 
is the finite component that we need to estimate. Approximating X^<lo w ' tn St^LoC/^ 5 )* ' s a wa y °^ ren dering 
the e-component finite while hopefully not altering the e x -component too much. There should be other substitutions 
that lead to better approximations (in this context, see the final paragraph in Section 1.1). 
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largest magnitude eigenvalue is 1 and has multiplicity 1, and the corresponding left eigenvector is 
the stationary distribution. We sort the eigenvalues A^ in decreasing order of magnitude, so that 
1 = Ai > | A2 1 > • • • > |A S | for some 2 < s < n. It turns out that | A2 1 determines the mixing time 
of the chain. 

The following theorem shows that if 1 — ft is small compared to 1 — | A2 1 , the gradient approx- 
imation described above is accurate. Since we will be using the estimate as a direction in which to 
update the parameters, the theorem compares the directions of the gradient and its estimate. In this 
theorem, k 2 {A) denotes the spectral condition number of a nonsingular matrix A, which is defined 
as the product of the spectral norms of the matrices A and A -1 , 

K 2 (A) = P|| 2 P _1 ||2, 

where 

IIAIU = max IIAeI], 

x:||x|| = l 

and 1 1 x ]| denotes the Euclidean norm of the vector x. 

Theorem 3. Suppose that the transition probability matrix P(9) satisfies Assumption 1 with sta- 
tionary distribution tt' = (m, ■ ■ ■ , n n ), and has n distinct eigenvalues. Let S = (xix? ■ ■ ■ x n ) be 
the matrix of right eigenvectors of P corresponding, in order, to the eigenvalues 1 = Ai > | A2 1 > 
• • • > Then the normalized inner product between V77 and fi'Vprj satisfies 

. - ^ s« liv^, v .^,ll ^_^_ m 

where II = diag(7Ti, . . . , n n ). 

Notice that rTlr is the expectation under the stationary distribution of r(X) 2 . 

As well as the mixing time (via \X2\X tne bound in the theorem depends on another parameter of 
the Markov chain: the spectral condition number of II 1 / 2 S. If the Markov chain is reversible (which 
implies that the eigenvectors xx, . . . ,x n are orthogonal), this is equal to the ratio of the maximum 
to the minimum probability of states under the stationary distribution. However, the eigenvectors 
do not need to be nearly orthogonal. In fact, the condition that the transition probability matrix 
have n distinct eigenvalues is not necessary; without it, the condition number is replaced by a more 
complicated expression involving spectral norms of matrices of the form (P — A^I). 

Proof. The existence of n distinct eigenvalues implies that P can be expressed as S'AS' -1 , where 
A = diag(Ai, . . . , A n ) (Lancaster & Tismenetsky, 1985, Theorem 4.10.2, p 153). It follows that for 
any polynomial /, we can write f(P) = Sf(A)S~ 1 . 

Now, Proposition 1 shows that Vrj — pVprj = V7r'(l — f$)Jp. But 

(1 - P)Jp = (1 - p) (r + (3Pr + (3 2 P 2 r + ■■■) 
= (1-P)(I + PP + @ 2 P 2 + ■ ■ ■ ) r 



P)sl J £p t A t )S- l r 



,t=o 



3=1 \t=o 
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where S 1 = {y x , y n )'. 

It is easy to verify that yi is the left eigenvector corresponding to A^, and that we can choose 
yi = 7r and x\ = e. Thus we can write 

n / oo \ 

(1 - 0)Jp = (1 - /3)en'r + £ ^ J> - fitfXtf r 



i=2 



< 1 -«-'+i>*i(i^) 



i=2 

(1 - /3)er? + SMS-\ 



where 



It follows from this and Proposition 1 that 



M = diag 0, 



1-/3 



1-/3 



V?7 ■ PV/sV 
||Vr/P 



1 



Vrj- (Vf7- Vtt'(1 -/9)Jg 



< 



HVi/P 

Vr; • Vtt' ((1 - ^)ery + SMS" x r) 

WW 

V'q ■ Vn'SMS^r 

WW 

IVtt'SMS-MI 



llVr/H 



by the Cauchy-Schwartz inequality. Since Vtt' = V y\fi?j II 1 / 2 , we can apply the Cauchy- 
Schwartz inequality again to obtain 



v 






|n 1 / 2 SM5~ 1 r| 




[Vull 



(28) 



We use spectral norms to bound the second factor in the numerator. It is clear from the definition 
that the spectral norm of a product of nonsingular matrices satisfies ||AB]|2 < ||v4]|2||-B]|2> an d that 
the spectral norm of a diagonal matrix is given by || diag(di, . . . ,rf n )]|2 = maxj \d{\. It follows that 



U 1/2 SMS 



n 1/2 5M5^ 1 n- 1/2 n 1/2 



< 



n l/2 5 S -l n -l/2 
2 



< k 2 (n 1/2 s) v^n^Y 

Combining with Equation (28) proves (27). 



2 

1-/8 



IMI 



/3|A 2 



□ 
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5. Estimating the Gradient in Parameterized Markov Chains 

Algorithm 1 introduces MCG (Markov Chain Gradient), an algorithm for estimating the approx- 
imate gradient Vpr) from a single on-line sample path Xq,X\, . . . from the Markov chain M(6). 
MCG requires only 2K reals to be stored, where K is the dimension of the parameter space: K 
parameters for the eligibility trace zt, and K parameters for the gradient estimate A t . Note that 
after T time steps At is the average so far of r(X t )z t , 

T-l 

Ar = -5>r(X t ). 

t=o 



Algorithm 1 The MCG (Markov Chain Gradient) algorithm 
l: Given: 

• Parameter 6 G M K . 

• Parameterized class of stochastic matrices V = {P(9) : 6 G M. K } satisfying Assumptions 
3 and 1. 

• /3e[o,i). 

• Arbitrary starting state Xg. 

• State sequence Xq,Xi, ... generated by M(9) (i.e. the Markov chain with transition 
probabilities P(6)). 

• Reward sequence r(Xo),r(X\), . . . satisfying Assumption 2. 

2: Set z = and A = (z , A G R K ). 
3: for each state Xt+i visited do 

4: z t+ i = 0Zt + ± j^- 

PX t X t+1 {8) 

5: A i+ i = A t + ^ [r{X t+l )z t+l - A t ] 
6: end for 



Theorem 4. Under Assumptions 1, 2 and 3, the MCG algorithm starting from any initial state Xq 
will generate a sequence Ao, Ai, . . . , At, . . . satisfying 

lim A t = Van w.p.l. (29) 

i-s>oo 

Proof. Let {X t } = {Xq,Xi, ...} denote the random process corresponding to M{6). If Xq ~ 7r 
then the entire process is stationary. The proof can easily be generalized to arbitrary initial distri- 
butions using the fact that under Assumption 1, {X t } is asymptotically stationary. When {X t } is 
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stationary, we can write 



>,3 



Y Pr(X 4 = i) Pv(X t+1 = j\X t = i)M E (j(i + l)\ Xt+l = j) 



1,3 



(30) 



where the first probability is with respect to the stationary distribution and J(t + 1) is the process 



j(t + i)= Y ^MXs). 

s=t+l 

The fact that E(J(i + = Jg(Xt+\) for all X t +\ follows from the boundedness of the 

magnitudes of the rewards (Assumption 2) and Lebesgue's dominated convergence theorem. We 
can rewrite Equation (30) as 



i,3 

where Xi(') denotes the indicator function for state i, 

Xi{Xt) 



Pij(0) 



J(t + l) 



1 ifX t = i, 
otherwise, 



and the expectation is again with respect to the stationary distribution. When X t is chosen according 
to the stationary distribution, the process {X t } is ergodic. Since the process {Z t } defined by 

Zt := Xi (X t ) Xj (X t+1 )^p^-J(t+l) 
Pij{V) 

is obtained by taking a fixed function of {X t }, {Z t } is also stationary and ergodic (Breiman, 1966, 



Proposition 6.31). Since 
(almost surely): 



is bounded by Assumption 3, from the ergodic theorem we have 



T-l 



tt'VP = V lim ^xdX t )x3(X t+ i)^j^J( 



'T-,oo T 



Pij{0) 



1 T_1 
lim — > 



t=0 
T-l 



lim — 



t=0 



Vpx t x t+1 (g) 
PXtXt+i (#) 

PXtXt+i (#) 



J(t + 1) 

T 00 

Y P M r{X s ) + Y P M r(Xs 

,s=i+l s=T+l 



(31) 
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Concentrating on the second term in the right-hand-side of (31), observe that: 

T-1 



1 E ^±^l E ^ r{Xs 



?U Px t x t+ M 



s=T+l 

T-1 

< 



1 T ~ X 



< 



T 

BR 



Vp Xt x t+1 (8) 



t=0 
T-1 oo 



PX t X t+1 (°) 



E ^-^Ir^, 

s=T+l 



E E z^ 1 



t=0 s=T+l 



si?^ (i - p T ) 

~ T(l -/3) 2 
-> Oas T -> oo, 

where R and S are the bounds on the magnitudes of the rewards and from Assumptions 2 

and 3. Hence, 

n'VPJ,= lim ly Vw ^ (fl) 



PX t X 1+1 (0) 



E /3 S ~*~V(X S ). 



(32) 



t=0 ' s=i+l 

Unrolling the equation for At in the MCG algorithm shows it is equal to 

T-1 _ , n s T 



' Px t x t+ M 
hence —> n'VPJp w.p.l as required. 



□ 



6. Estimating the Gradient in Partially Observable Markov Decision Processes 

Algorithm 1 applies to any parameterized class of stochastic matrices P(0) for which we can com- 
pute the gradients Vpij{0). In this section we consider the special case of P{6) that arise from a 
parameterized class of randomized policies controlling a partially observable Markov decision pro- 
cess (POMDP). The 'partially observable' qualification means we assume that these policies have 
access to an observation process that depends on the state, but in general they may not see the state. 

Specifically, assume that there are N controls U = {l,...,iV} and M observations y = 
{1, . . . , M}. Each u G U determines a stochastic matrix P{u) which does not depend on the 
parameters 6. For each state i G S, an observation Y G y is generated independently according to 
a probability distribution v(i) over observations in y. We denote the probability of observation y 
by v y {i). A randomized policy is simply a function \i mapping observations y G y into probability 
distributions over the controls U. That is, for each observation y, (i(y) is a distribution over the 
controls in U. Denote the probability under [x of control u given observation y by /j, u {y)- 

To each randomized policy //(■) and observation distribution u(-) there corresponds a Markov 
chain in which state transitions are generated by first selecting an observation y in state % according 
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to the distribution then selecting a control u according to the distribution fi(y), and then gen- 
erating a transition to state j according to the probability Pij(u). To parameterize these chains we 
parameterize the policies, so that p, now becomes a function fj,(6, y) of a set of parameters 6 6 M. K as 
well as the observation y. The Markov chain corresponding to 6 has state transition matrix \pij(8)] 
given by 

Pij(6) = 'E Y ^(i)'Eur~ JfJ ,(0,Y)Pij (U). (33) 

Equation (33) implies 

VpiiW = Vy(})pij(u)Vfi u (0, y). (34) 

Algorithm 2 introduces the GPOMDP algorithm (for Gradient of a Partially Observable Markov 
Decision Process), a modified form of Algorithm 1 in which updates of zt are based on ^Uti^i Yt)> 
rather than px t x t+1 {6)- Note that Algorithm 2 does not require knowledge of the transition prob- 
ability matrix P, nor of the observation process u; it only requires knowledge of the randomized 
policy /j,. GPOMDP is essentially the algorithm proposed by Kimura et al. (1997) without the 
reward baseline. 

The algorithm GPOMDP assumes that the policy \x is a function only of the current observation. 
It is immediate that the same algorithm works for any finite history of observations. In general, an 
optimal policy needs to be a function of the entire observation history. GPOMDP can be extended 
to apply to policies with internal state (Aberdeen & Baxter, 2001). 

Algorithm 2 The GPOMDP algorithm, 
l: Given: 

• Parameterized class of randomized policies {/j,(6, ■) : 9 6 K } satisfying Assumption 4. 

• Partially observable Markov decision process which when controlled by the randomized 
policies /u(#, -) corresponds to a parameterized class of Markov chains satisfying As- 
sumption 1. 

• £e[0,l). 

• Arbitrary (unknown) starting state Xq. 

• Observation sequence Yq,Yi, . . . generated by the POMDP with controls Uq, U\, . . . 
generated randomly according to /j,(8,Y t ). 

• Reward sequence r(Xo), r(X\), . . . satisfying Assumption 2, where Xo, X\, ... is the 
(hidden) sequence of states of the Markov decision process. 

2: Set z = and A = (z , A G R K ). 

3: for each observation Y t , control Ut, and subsequent reward r(Xt+i) do 
4: z t+ i = pzt + , 

Vu t (8,Y t ) 

5: A m = A t + [r(X t+1 )z t+1 - A t ] 
6: end for 
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For convergence of Algorithm 2 we need to replace Assumption 3 with a similar bound on the 
gradient of \x: 



Assumption 4. The derivatives, 



dHu{0,y) 
dBk 



exist for all u £U, y G y and 6 6 M . The ratios 



y=l...M;u=l...N;k=l...K 



are uniformly bounded by < oo for all 6 6 M . 

Theorem 5. Under Assumptions 1, 2 and 4, Algorithm 2 starting from any initial state Xq will 
generate a sequence Aq, Ai, . . . , Af, . . . satisfying 



lim At = Vrt] w.p.l. 

t— ¥00 

Proof. The proof follows the same lines as the proof of Theorem 4. In this case, 

= ^2 Tr(i)pij(u)vy(i)Vn u (6,y)Jp(j) from (34) 
V/i«(e s y) 



(35) 



■iiu(0,y)Jp(j), 



where the expectation is with respect to the stationary distribution of {X t }, and the process {Z' t } is 
defined by 

Z' t := Xi (X t ) Xj (X t+1 ) Xn (U t ) X y(Y t ) Vt " u f' f j(t + l), 

Hu\v-,y) 

where Ut is the control process and Y t is the observation process. The result follows from the same 
arguments used in the proof of Theorem 4. □ 

6.1 Control dependent rewards 

There are many circumstances in which the rewards may themselves depend on the controls u. For 
example, some controls may consume more energy than others and so we may wish to add a penalty 
term to the reward function in order to conserve energy. The simplest way to deal with this is to 
define for each state i the expected reward r(i) by 



r(i) = Ey^E^^r ([/,?), 



(36) 
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and then redefine J« in terms of r: 



Jf}(8,i) := lim E 6 

N— >oo 



N 



t=0 



An 



(37) 



where the expectation is over all trajectories Xq , X\ , . . . . The performance gradient then becomes 

\?n = Vyr'f + yr'Vf, 

which can be approximated by 

Vn'n = vr' [vpJ/3 + Vf] , 

due to the fact that Jp satisfies the Bellman equations (20) with f replaced by r. 

For GPOMDP to take account of the dependence of r on the controls, its fifth line should be 
replaced by 

1 



A 



t+i 



t + 1 



r{U t +i,X t +i) zt+i + — 



A, 



It is straightforward to extend the proofs of Theorems 2, 3 and 5 to this setting. 
6.2 Parameter dependent rewards 

It is possible to modify GPOMDP when the rewards themselves depend directly on 9. In this case, 
the fifth fine of GPOMDP is replaced with 

1 

At+i 



A* + 



t + 1 



[r(6,X t+l )z t+1 + Vr(6,X t+l ) - A t 



(38) 



Again, the convergence and approximation theorems will carry through, provided Vr(0, i) is uni- 
formly bounded. Parameter-dependent rewards have been considered by Glynn (1990), Marbach 
and Tsitsiklis (1998), and Baird and Moore (1999). In particular, Baird and Moore (1999) showed 
how suitable choices of r(6,i) lead to a combination of value and policy search, or "VAPS". For 
example, if J(6, i) is an approximate value-function, then setting 13 



r{6,X t ,X t ^) 



r(X t )+aJ(e,X t )-J(e,X t ^) 



where r(X t ) is the usual reward and a 6 [0, 1) is a discount factor, gives an update that seeks to 
minimize the expected Bellman error 

2 



E 7 ^ 



i=l 



J(9,i) 



(39) 



This will have the effect of both minimizing the Bellman error in J(6, i), and driving the system 
(via the policy) to states with small Bellman error. The motivation behind such an approach can 
be understood if one considers a J that has zero Bellman error for all states. In that case a greedy 
policy derived from J will be optimal, and regardless of how the actual policy is parameterized, the 
expectation of Ztr{6, Xt, Xt-i) will be zero and so will be the gradient computed by GPOMDP. 
This kind of update is known as an actor-critic algorithm (Barto et al, 1983), with the policy playing 
the role of the actor, and the value function playing the role of the critic. 



13. The use of rewards r(B,Xt,Xt-i) 
analysis. 



that depend on the current and previous state does not substantially alter the 
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6.3 Extensions to infinite state, observation, and control spaces 

The convergence proof for Algorithm 2 relied on finite state (S), observation (30 and control (U) 
spaces. However, it should be clear that with no modification Algorithm 2 can be applied imme- 
diately to POMDPs with countably or uncountably infinite S and y, and countable U. All that 
changes is that Pij(u) becomes a kernel p(x, x', u) and u(i) becomes a density on observations. In 
addition, with the appropriate interpretation of V/u/^i, it can be applied to uncountable U. Specifi- 
cally, if U is a subset of M. N then /j,(y, 6) will be a probability density function on U with fj, u (y, 0) 
the density at u. If IA and y are subsets of Euclidean space (but S is a finite set), Theorem 5 can be 
extended to show that the estimates produced by this algorithm converge almost surely to \/gr/. In 
fact, we can prove a more general result that implies both this case of densities on subsets of M N as 
well as the finite case of Theorem 5. We allow U and y to be general spaces satisfying the following 
topological assumption. (For definitions see, for example, (Dudley, 1989).) 

Assumption 5. The control space U has an associated topology that is separable, Hausdorff, and 
first-countable. For the corresponding Borel o-algebra B generated by this topology, there is a 
o -finite measure A defined on the measurable space (U,B). We say that A is the reference measure 
forU. 

Similarly, the observation space y has a topology, Borel o-algebra, and reference measure 
satisfying the same conditions. 

In the case of Theorem 5, where U and y are finite, the associated reference measure is the 
counting measure. For U = M N and y = M. M , the reference measure is Lebesgue measure. We 
assume that the distributions v(i) and p,(6, y) are absolutely continuous with respect to the reference 
measures, and the corresponding Radon-Nikodym derivatives (probability masses in the finite case, 
densities in the Euclidean case) satisfy the following assumption. 

Assumption 6. For every y E y and 8 E R , the probability measure [i(8, y) is absolutely contin- 
uous with respect to the reference measure for U. For every i E S, the probability measure v(i) is 
absolutely continuous with respect to the reference measure for y. 

Let A be the reference measure for U. For all u E U, y E y, 6 E R , and k E {1, . . . , K}, the 
derivatives 

d dfi(6,y) 
d6 k dX [u) 

exist and the ratios 

a dfi n (0,y) / \ 



dnu(0,y) 

are bounded by < oo. 



With these assumptions, we can replace \i in Algorithm 2 with the Radon-Nikodym derivative 
of fj, with respect to the reference measure on U. In this case, we have the following convergence 
result. This generalizes Theorem 5, and also applies to densities ^ona Euclidean space U. 

Theorem 6. Suppose the control space U and the observation space y satisfy Assumption 5 and let 
A be the reference measure on the control space U. Consider Algorithm 2 with 
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replaced by 



Under Assumptions 1, 2 and 6, this algorithm, starting from any initial state Xq will generate a 
sequence Ao, Ai, . . . , At, . . . satisfying 

lim A t = Vrtj w.p.l. 
t— too 

Proof. See Appendix B □ 



7. New Results 

Since the first version of this paper, we have extended GPOMDP to several new settings, and also 
proved some new properties of the algorithm. In this section we briefly outline these results. 



7.1 Multiple Agents 

Instead of a single agent generating actions according to /z(0, y), suppose we have multiple agents 
% = 1, . . . , n a , each with their own parameter set 6 l and distinct observation of the environment 
y\ and that generate their own actions u l according to a policy /j, u i(9 t ,y l ). If the agents all re- 
ceive the same reward signal r(X t ) (they may be cooperating to solve the same task, for example), 
then GPOMDP can be applied to the collective POMDP obtained by concatenating the observa- 
tions, controls, and parameters into single vectors y = [y 1 , . . . , y n "] , u = [u 1 , . . . , it" a ] , and 
6 = [9 1 , . . . , 8 na ] respectively. An easy calculation shows that the gradient estimate A generated 
by GPOMDP in the collective case is precisely the same as that obtained by applying GPOMDP to 
each agent independently, and then concatenating the results. That is, A = [A 1 , . . . , A" a ] , where 
A 1 is the estimate produced by GPOMDP applied to agent i. This leads to an on-line algorithm 
in which the agents adjust their parameters independently and without any explicit communication, 
yet collectively the adjustments are maximizing the global average reward. For similar observa- 
tions in the context of REINFORCE and VAPS, see Peshkin et al. (2000). This algorithm gives a 
biologically plausible synaptic weight-update rule when applied to networks of spiking neurons in 
which the neurons are regarded as independent agents (Bartlett & Baxter, 1999), and has shown 
some promise in a network routing application (Tao, Baxter, & Weaver, 2001). 



7.2 Policies with internal states 

So far we have only considered purely reactive or memoryless policies in which the chosen control 
is a function of only the current observation. GPOMDP is easily extended to cover the case of 
policies that depend on finite histories of observations Y tl it_i, . . . , Y t -k, but in general, for optimal 
control of POMDPs, the policy must be a function of the entire observation history. Fortunately, the 
observation history may be summarized in the form of a belief state (the current distribution over 
states), which is itself updated based only upon the current observation, and knowledge of which 
is sufficient for optimal behaviour (Smallwood & Sondik, 1973; Sondik, 1978). An extension of 
GPOMDP to policies with parameterized internal belief states is described by Aberdeen and Baxter 
(2001), similar in spirit to the extension of VAPS and REINFORCE described by Meuleau et al. 
(1999). 
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7.3 Higher-Order Derivatives 

GPOMDP can be generalized to compute estimates of second and higher-order derivatives of the 
average reward (assuming they exist), still from a single sample path of the underlying POMDP. 
To see this for second-order derivatives, observe that if rj(9) = f q(9,x)r(x) dx for some twice- 
differentiable density q(0, x) and performance measure r(x), then 

2, 



V 2 ^) = [r(x)^^± q (6,x)dx 



where V 2 denotes the matrix of second derivatives (Hessian). It can be verified that 

V 2 q{9,x) 



q(9,x) 



V 2 \ogq(6,x) + [V\ogq{6,x)Y (40) 



where the second term on the right-hand-side is the outer product between V log q(6, x) and itself 
(that is, the matrix with entries d/dOi log q(8, x)d/d9j log q(6, x)). Taking a; to be a sequence of 
states Xq, X\, . . . , Xt between visits to a recurrent state i* in a parameterized Markov chain (recall 
Section 1.1.1), we have q(9, X) = HjSQPx t x t+ i ($)> which combined with (40) yields 



q(6,X) 



T-1 

E 

t=0 



V 2 Px t x t+1 (9) 
Px t x t+1 (0) 



T-1 

E 



Px t x t+1 (0) 



+ 



T-1 

E 

.t=o 



Vpx t x t+1 (g) 
Px t x t+1 (0) 



-i £ 



(the squared terms in this expression are also outer products). From this expression we can derive 
a GPOMDP-like algorithm for computing a biased estimate of the Hessian V 2 r/(#), which involves 
maintaining — in addition to the usual eligibility trace z t — a second matrix trace updated as follows: 



H+l 



V 2 Px t x t+1 {9) 
Px t x t+1 (0) 



Vpx t x t+1 (9) 
Px t x t+1 ( e ) 



After T time steps the algorithm returns the average so far of r(X t ) [Z t + zf] where the second term 
is again an outer product. Computation of higher-order derivatives could be used in second-order 
gradient methods for optimization of policy parameters. 



7.4 Bias and Variance Bounds 

Theorem 3 provides abound on the bias of Vpr]{9) relative to Vrj(6) that applies when the underly- 
ing Markov chain has distinct eigenvalues. We have extended this result to arbitrary Markov chains 
(Bartlett & Baxter, 2001). However, the extra generality comes at a price, since the latter bound in- 
volves the number of states in the chain, whereas Theorem 3 does not. The same paper also supplies 
a proof that the variance of GPOMDP scales as 1/(1 — /3) 2 , providing a formal justification for the 
interpretation of j3 in terms of bias/variance trade-off. 



8. Conclusion 

We have presented a general algorithm (MCG) for computing arbitrarily accurate approximations 
to the gradient of the average reward in a parameterized Markov chain. When the chain's transition 
matrix has distinct eigenvalues, the accuracy of the approximation was shown to be controlled by the 
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size of the subdominant eigenvalue [ A2 ] . We showed how the algorithm could be modified to apply 
to partially observable Markov decision processes controlled by parameterized stochastic policies, 
with both discrete and continuous control, observation and state spaces (GPOMDP). For the finite 
state case, we proved convergence with probability 1 of both algorithms. 

We briefly described extensions to multi-agent problems, policies with internal state, estimating 
higher-order derivatives, generalizations of the bias result to chains with non-distinct eigenvalues, 
and a new variance result. There are many avenues for further research. Continuous time results 
should follow as extensions of the results presented here. The MCG and GPOMDP algorithms can 
be applied to countably or uncountably infinite state spaces; convergence results are also needed in 
these cases. 

In the companion paper (Baxter et al., 2001), we present experimental results showing rapid 
convergence of the estimates generated by GPOMDP to the true gradient Vr?. We give on-line 
variants of the algorithms of the present paper, and also valiants of gradient ascent that make use of 
the estimates of Vp'q. We present experimental results showing the effectiveness of these algorithms 
in a variety of problems, including a three-state MDP, a nonlinear physical control problem, and a 
call-admission problem. 



Acknowledgements 

This work was supported by the Australian Research Council, and benefited from the comments of 
several anonymous referees. Most of this research was performed while the authors were with the 
Research School of Information Sciences and Engineering, Australian National University. 



Appendix A. A Simple Example of Policy Degradation in Value-Function Learning 



Approximate value-function approaches to reinforcement work by minimizing some form of error 
between the approximate value function and the true value function. It has long been known that this 
may not necessarily lead to improved policy performance from the new value function. We include 
this appendix because it illustrates that this phenomenon can occur in the simplest possible system, 
a two-state MDP, and also provides some geometric intuition for why the phenomenon arises. 

Consider the two-state Markov decision process (MDP) in Figure 1. There are two controls 
«i, U2 with corresponding transition probability matrices 



P{u x ) 



P{u 2 ) 



so that u\ always takes the system to state 2 with probability 2/3, regardless of the starting state (and 
therefore to state 1 with probability 1/3), and u 2 does the opposite. Since state 2 has a reward of 1, 
while state 1 has a reward of 0, the optimal policy is to always select action u\. Under this policy 
the stationary distribution on states is [7^,7^] = [1/3,2/3], while the infinite-horizon discounted 
value of each state i = 1, 2 with discount value a G [0, 1) is 



Ja(i) 




c?r{X t ) 



where the expectation is over all state sequences Xq , X\ , X 2 , ■ ■ ■ with state transitions generated ac- 
cording to P(u\). Solving Bellman's equations: J a = r + aP(u\)J a , where J a = [J Q (1), J Q (2)]' 
and r = [r(l),r(2)] ; yields J a (l) = 3^ and J a (2) = 1 + 



2a 



3(l-a) ' 
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r(1) = 



r(2) = 1 




Figure 1 : Two-state Markov Decsision Process 



Now, suppose we are trying to learn an approximate value function J for this MDP, i.e. , J(i) = 
W(f>(i) for each state i = 1, 2 and some scalar feature (f) must have dimensionality 1 to ensure that 
J really is approximate). Here w G K is the parameter being learnt. For the greedy policy obtained 
from J to be optimal, J must value state 2 above state 1. For the puiposes of this illustration choose 
(f)(1) = 2, cj){2) = 1, so that for j(2) > J(l), w must be negative. 

Temporal Difference learning (or TD(A)) is one of the most popular techniques for training 
approximate value functions (Sutton & Barto, 1998). It has been shown that for linear functions, 
TD(1) converges to a parameter w* minimizing the expected squared loss under the stationary 
distribution (Tsitsikilis & Van-Roy, 1997): 



w 



argmin^ ) tt; [w<p(i) - J a (i)Y 



i=l 



(41) 



Substituting the previous expressions for 7ri,7T2,</> and J a under the optimal policy and solving 
for w*, yields w* = gnz?t) 1 Hence w* > for all values of a G [0, 1), which is the wrong 
sign. So we have a situation where the optimal policy is implementable as a greedy policy based 
on an approximate value function in the class (just choose any w < 0), yet TD(1) observing the 
optimal policy will converge to a value function whose corresponding greedy policy implements the 
suboptimal policy. 

A geometrical illustration of why this occurs is shown in Figure 2. In this figure, points on the 
graph represent the values of the states. The scales of the state 1 and state 2 axes are weighted by 
and \/n(2) respectively. In this way, the squared euclidean distance on the graph between 
two points J and J corresponds to the expectation under the stationary distribution of the squared 
difference between values: 



V^(T)J(1),V^(2)J(2) - J(l), V^(2) J(2) 



E, [J(X)-J(X) 



For any value function in the shaded region, the corresponding greedy policy is optimal, since 
those value functions rank state 2 above state 1. The bold line represents the set of all realizable 
approximate value functions (w (f)(1), w (f)(2)). The solution to (41) is then the approximate value 
function found by projecting the point corresponding to the true value function [(J a (l), J a (2)] onto 
this line. This is illustrated in the figure for a = 3/5. The projection is suboptimal because weighted 
mean-squared distance in value-function space does not take account of the policy boundary. 



Appendix B. Proof of Theorem 6 

The proof needs the following topological lemma. For definitions see, for example, (Dudley, 1989, 
pp. 24-25). 
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Figure 2: Plot of value-function space for the two-state system. Note that the scale of each axis has 
been weighted by the square root of the stationary probability of the corresponding state 
under the optimal policy. The solution found by TD(1) is simply the projection of the true 
value function onto the set of approximate value functions. 



Lemma 7. Let (X, T) be a topological space that is Hausdorjf, separable, and first-countable. 
Let B be the Borel a-algebra generated by T. Then the measurable space (X, B) has a sequence 
Si, <?2j ■ ■ ■ B of sets that satisfies the following conditions: 

1. Each Si is a partition of X (that is, X = y}{S: S G Si} and any two distinct elements of Si 
have empty intersection). 

2. For all x G X, {x} G B and 

oo 

f]{S e S t : x e S} = {x}. 

Proof. Since X is separable, it has a countable dense subset S = {x\,X2, . . .}. Since X is first- 
countable, each of these Xi has a countable neighbourhood base, iVj. Now, construct the partitions 
Si using the countable set N = (J^i ^ as follows. Let So = X and, for i = 1, 2, . . ., define 

Si = {S n Ni : S G U{Sfl (X — Ni) : S G . 
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Clearly, each Si is a measurable partition of X. Since X is Hausdorff, for each pair x, x! of distinct 
points from X, there is a pair of disjoint open sets A and A' such that x £ A and x' £ yl'. Since S 1 
is dense, there is a pair s, 5' from S with 5 £ A and s' £ ^4'. Also, N contains neighbourhoods N s 
and N s i with iV s C A and iV s / C ^4'. So N s and jV s / are disjoint. Thus, for sufficiently large i, x 
and x' fall in distinct elements of the partition S{. Since this is true for any pair x, x', it follows that 

00 

p|{5 £ S t : x £ S} C {a;}. 

i=l 

The reverse inclusion is trivial. The measurability of all singletons {x} follows from the measura- 
bility of S x := [j^S £ S t : S n = </>} and the fact that {x} = X - S x . □ 

We shall use Lemma 7 together with the following result to show that we can approximate 
expectations of certain random variables using a single sample path of the Markov chain. 

Lemma 8. Let (X, B) be a measurable space satisfying the conditions of Lemma 7, and let S\ , S2 , • • • 
be a suitable sequence of partitions as in that lemma. Let fibe a probability measure defined on this 
space. Let f be an absolutely integrable function on X. For an event S, define 

HS) = i^. 

For each x £ X and k = 1,2,..., let Sk(x) be the unique element of S^ containing x. Then for 
almost all x in X, 

lim f(S k (x)) = /Or). 

k — loo 

Proof. Clearly, the signed finite measure cf) defined by 

</>{E) = [ fdfjt (42) 
Je 

is absolutely continuous with respect to \i, and Equation (42) defines / as the Radon-Nikodym 
derivative of <j) with respect to \x. This derivative can also be defined as 

**(x) = lim 

d/i fc-nxi fl(Sk(x)) 

See, for example, (Shilov & Gurevich, 1966, Section 10.2). By the Radon-Nikodym Theorem (Dud- 
ley, 1989, Theorem 5.5.4, p. 134), these two expressions are equal a.e. (^). □ 

Proof. (Theorem 6.) From the definitions, 

Vp-q = tt'VPJ/3 

n n 

.=1 i=i 

For every y, \x is absolutely continuous with respect to the reference measure A, hence for any % and 
j we can write 

dn{Q,y) 



Pij{6)= / / Pij{u) — tt — {u)d\{u)dv{i){y). 



iyJu d\ 
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Since A and v do not depend on and dp(6,y)/d\ is absolutely integrable, we can differentiate 
under the integral to obtain 



Vpy(0) =11 KjWV 
Jy Ju 



dX 



u) d\(u) du(i)(y). 



To avoid cluttering the notation, we shall use fj, to denote the distribution fi(6,y) on U, and v to 
denote the distribution v{%) on y. With this notation, we have 



r r V— 

Vpij(^) = / / Pij^dndu 
JyJu £ 



Now, let p be the probability measure onJ'xW generated by jj, and za We can write (43) as 



V ^ = 22 W) / Pij^r d P- 



Using the notation of Lemma 8, we define 



Pij(S) 



V(S) 



fsPij d P 
P(S) ' 

P(S) Is f x 



for a measurable set S C ^ x jy/. Notice that, for a given ?, j, and S 1 , 



Pii(^) = Pr (^t+i = 3 \ x t = h {y,u) 6 S) 



( 

V(5) = E ' dX 



du 



X t = i, (Y t , C/*) 6 S . 



Let S\, <?2) ■ ■ ■ be a sequence of partitions of 3^ x W as in Lemma 7, and let <Si,(y, u) denote the 
element of Sk containing (y, u). Using Lemma 8, we have 



r \7ilL 
/ v — ^ 



lim Pij(S k (y,u)) V (S k {y,u)) dp(y,u) 



k— >oo 



lim £ / Ki (5)V(5)dp, 



ses k 
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where we have used Assumption 6 and the Lebesgue dominated convergence theorem to interchange 
the integral and the limit. Hence, 



= lim Y, E ^)p(S)Pv(S)Mj)V(S) 
k— too z — ' f — ; 



y ses k 



lim V Pi(X t = i) Pr((Y t , C/ t ) G 5) Pr = j \X t = i, (Y t , U t ) G S) 

k—xx> z — ' 



V(J(t+l)\X t+l =j)E 



lim VE 



i,j,S 



X t = i, (Y u U t ) G S 



Xi{X t )xs{YuU t ) X j{X t+l )J{t + l)- 



dx 



dX 



where probabilities and expectations are with respect to the stationary distribution it of X t , and the 
distributions on Y t , Ut- Now, the random process inside the expectation is asymptotically stationary 
and ergodic. From the ergodic theorem, we have (almost surely) 



V 0V = lim lim 1 J2 E Xi(X t )xs(Y t , U t ) Xj (X t+1 )J(t + 1) 

k— ¥oo T— too 1 *■ — • *■ — ' 

i,j,S t=0 

It is easy to see that the double limit also exists when the order is reversed, so 

T-l 



T-l 



V 



dfi 
dX 



d/j. 
dX 



i * 1 V^H 

V PV = lim - E i im E^( x ')^(y,^)x.(^m)^ + l)^f 

T— foo J z — ' k — yoc z — ' iiii 
*=0 i,f,S dX 



t=0 i,j,S 

T-l yjdjjJ^Yt) 



ly V^r(Ut) 
T ^ du(6,Y t ) itt \ v y 



t=0 dX \ u t) 

The same argument as in the proof of Theorem 4 shows that the tails of J(t + 1) can be ignored 
when 



, dfx{e,Y t ) 
dX 



(Ut 



and |r"(Xf)| are uniformly bounded. It follows that — > n'WPJp w.p. 1, as required. 



□ 
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